Polarization and localization in insulators: Generating function approach 



Ivo Souza, Tim Wilkens, and Richard M. Martin 
Department of Physics and Materials Research Laboratory, Umversity of Ilhnois, Urbana, Illinois 61801 

(February 1, 2008) 

We develop the theory and practical expressions for the full quantum-mechanical distribution of 
the intrinsic macroscopic polarization of an insulator in terms of the ground state wavefunction. 
The central quantity is a cumulant generating function which yields, upon successive differentiation, 
all the cumulants and moments of the probability distribution of the center of mass 'X./N of the 
electrons, defined appropriately to remain valid for extended systems obeying twisted boundary con- 
ditions. The first moment is the average polarization, where we recover the well-known Berry phase 
expression. The second cumulant gives the mean-square fiuctuation of the polarization, which defines 
an electronic localization length along each direction i: = {{Xf) — {Xi)'^)/N. It follows from 
the fluctuation-dissipation theorem that in the thermodynamic limit diverges for metals and is a 
flnite, measurable quantity for insulators. It is possible to define for insulators maximally-localized 
"many-body Wannier functions" , which for large A'^ become localized in disconnected regions of the 
high-dimensional configuration space, establishing a direct connection with Kohn's theory of the 
insulating state. Interestingly, the expression for , which involves the second derivative of the 
wavefunction with respect to the boundary conditions, is directly analogous to Kohn's formula for 
the "Drude weight" as the second derivative of the energy. 



An insulator is distinguished from a conductor at zero temperature by its vanishing dc conductivity and its ability 
to sustain a macroscopic polarization, both with and without an applied electric field In the classical theory of 

electromagnetism in materials this distinction is often cast in terms of the difference between "free charges" that carry 
the dc current in a conductor and polarizable "bound charges" in an insulator |^,^ . Such a description conflicts with 
the fact that, even in highly ionic solids, the electrons are not well localized near the ions, and there is appreciable 
interpenetration between the ionic charge densities Q-Q] ■ The inadequacy of such a textbook picture is particularly 
striking in the case of covalent insulators, whose charge density is delocalized, as in metals Therefore the qualitative 
difference between metals and insulators is not apparent from inspection of the charge distribution, and the correct 
notion of electronic localization in insulators versus delocalization in conductors must be sought elsewhere. 

As shown by Kohn [^|J^, localization is a property of the many-electron wavefunction: insulating behavior arises 
whenever the ground state wavefunction of an extended system breaks up into a sum of functions which are 
localized in essentially disconnected regions TZm. of the high-dimensional configuration space. When using periodic 
boundary conditions on a supercell containing N electrons, an insulating wavefunction can be written as 



where for a large supercell and vI/m' have an exponentially small overlap for M' ^ M. Hence, electronic 
localization in insulators does not occur in real space (charge density) but in configuration space (wavefunction). 
Kohn argued that such disconnectedness is in fact the signature of an insulating wavefunction. 

The other aspect of insulators, macroscopic polarization, is a subject about which there has been much debate. In a 
finite system with open boundary conditions, it is simply the dipole moment of the charge distribution divided by the 
volume In the case of an extended system, such as a crystalline solid, the situation is far less clear: although there 
exist well-known expressions for dielectric response functions [ p^ , pl| , the very definition of macroscopic polarization 
as a bulk property independent of surface termination remained controversial for a long time. Only in recent years has 
a theory of polarization emerged - the so-called Berry phase formulation |p"2|-|r7[ - for the average bulk polarization 
in terms of the ground state wavefunction of an insulating crystal. This theory shows that in general the information 
about the macroscopic polarization of a periodic system is not in the charge density, but in the wavefunction. This 
important finding is consistent with the well-known fact that the dipole moment of a periodic, continuous charge 
distribution is ill-defined since the expressions for the first moment of the charge distribution are valid for 

finite systems but do not have a well-defined thermodynamic limit independent of the surface. The only instance 
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where a meaningful dipole moment can be assigned to a unit cell of the crystal is when the charge distribution in 
the unit cell can be resolved into contributions which are localized in non-overlapping regions and can be ascribed 
to identifiable atoms (ions, molecules) (the so-called "Clausius-Mossotti limit"). However, such a limit is rather 
unrealistic for most insulators, with the possible exception of some organic and molecular crystals. 

In this paper we present a comprehensive theory of polarization and localization in insulators, which generalizes 
the Berry phase theory and merges with Kohn's theory of localization in the insulating state. By introducing the 
generating function formalism, we show the rigorous connections between polarization and localization, establish 
relations to experimentally measureable quantities, and provide formulas for practical calculations. 

The paper is organized as follows: in Section || we summarize some of the recent developments in the field, pointing 
to some of the open issues which will be addressed in this work, and collect relevant equations for later reference. 
In Section III is introduced the central concept in the present formulation: the generating function applied to the 
quantum probability distribution of the macroscopic polarization. The simple case of a finite system is discussed first, 
and then we show how to generalize the expressions to deal with extended systems. In Section |^ we establish the 
fluctuation-dissipation relation between the quadratic quantum fluctuations of the polarization and the absorptive 
part of the conductivity for an extended system. In Section we discuss the qualitative differences between the 
quadratic fluctuations in insulators and in conductors in terms of a properly defined localization length; for insulators 
we derive an inequality involving the localization length and the optical gap. In Section VI we introduce "many-body 
Wannier functions" and establish the connection between the present formalism and Kohn's theory of localization 
in the insulating state. Discretized formulas for the polarization and localization which can be used in numerical 



many-body calculations are derived in Section VIl . A discussion of our results is presented in Section VIII 



II. OVERVIEW OF RECENT DEVELOPMENTS 



A meaningful definition of average macroscopic polarization which is generally applicable - to both finite and 
extended systems - and relates closely to the way polarization is experimentally measured, can be obtained by taking 
as the more basic concept the change in polarization induced by a slow change in some parameter A in the hamiltonian 

|^,|ll. The resulting expressions for an extended system are in terms of the derivative dF^^J^^/dX, and Resta ]2l| ] 
proposed to calculate the finite change in bulk polarization as 



APbuik = / dX^^. (2) 
Jo 9X 

Considering the variation in A as an adiabatic time evolution, dP\^Ji^/ dX is the spatially-averaged adiabatic current 
flowing through the bulk; thus this equation gives APbuik as an integrated bulk current. According to classical 
electrodynamics [H,0 , this adiabatic polarization current is 



ap[,l 1 r . .a. 1 



rfrja-T^J^'\ (3) 



dx y./v V 

where V is the volume of the system and Jbulk current density in the bulk, where A plays the role of time in the 

usual expression for the current (for a derivation see for example Ref. jl^). J*-'^-' can be expressed as the adiabatic 
limit of a Kubo formula for the current [^l|j2^ ]. In its usual form the Kubo formula involves a summation over all 
eigenstates (ground state and excited states), which not only makes it impractical for actual calculations, but also 
fails to reflect the fact that bulk polarization is a ground state property. 

In the Berry phase theory of polarization p^-p^, the Kubo formula is recast in a form which only depends on 
the ground state wavefunction; in order to arrive at such an expression, it is convenient to impose twisted boundary 
conditions over the volume V on the many-body wavefunction, which becomes labelled by k: 



<^'(xi, ...,x, + L, ...,XAr) = e'''-L*^^\xi, ...,x„ ...,x^), (4) 

where N is the number of electrons in the system, L ~ (Li,0, 0), for example, and — 7r/ij < ki < n/Li. It is 
convenient to introduce the following wavefunction: 
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where X = X^^j^ x^, so that 'K/N is the position operator for the center of mass of the N electrons in the volume V. 

^[^^ can be regarded as the many-body analogue of the cell-periodic part of the Bloch function in the single-electron 
case. It obeys periodic boundary conditions, and the k-dependence is transfered from the boundary conditions to the 
hamiltonian: if "^[^^ is the ground state of an hamiltonian H^'^\ then is the ground state of the hamiltonian 



(A)g»k-X 



(6) 



which, for a non-relativistic hamiltonian without spin-orbit coupling, can be obtained from H^^^ by performing the 
gauge transformation Pi ^ Pi -|- hh on the momentum operator of each particle (for the single-particle analogue, see 
Ref. [^.) Using the function the Kubo formula for the adiabatic current J^^-* for a particular choice of twisted 
boundary conditions can be expressed as p4,E3,E6| 



(7) 



which indeed only depends on the ground state wavefunction. Substituting this expression into Eq. |^ and using Eq. 
0, it can be shown that the net change in polarization along the path parametrized by A is 



(APel 
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dk I d\ 
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where the integral in k is over all twisted boundary conditions. 

As expected for a measurable quantity, the above expression, as well as Eq. |^, are invariant under gauge transfor- 
mations of the form 



(9) 



where (p{k, A) is a smooth, real function. The ground state wavefunctions at k and k + G, where G is a basis 
vector of the reciprocal lattice of the cell of volume V, can differ at most by a global phase factor, since the boundary 
conditions are the same and the ground state is assumed to be non-degenerate |p5[ : 



(A) ^ pie(k,A;G)^j,_(A) 



k+G- 
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If 0(k, A; G), which is at our disposal, is chosen to be independent of A, then it can be shown that the net change in 
polarization becomes simply 



p(0) 
-^el ' 



(11) 



where 



,(A) 
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(12) 



and again the integral in k is over all twisted boundary conditions. Eqs. |[ ^ and ^ form the central result of 
Ref. |14| , which gives the many-body generalization of the Berry phase theory of polarization, originally formulated 
by King-Smith and Vanderbilt for independent electrons p^ , p^ . In the derivation it is assumed that the ground 
state is isolated from the excited states by a finite energy gap, and that there are no long-range correlations [ p5| . 
In the independent-electron theory the Berry phase formula in terms of Bloch functions can be recast in terms of 
localized Wannier functions 129], yielding the intuitive result that the electronic polarization is given by the sum of 
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the centers of charge of the occupied Wa nnier functions ||T2|,|l3[. Unfortunately, Wannier functions are only defined 
in a one-electron framework. In Section VI B we will introduce functions which are in some sense the many-body 
counterpart of Wannier functions (we will term them "many-body Wannier functions"), and can be chosen to be 
localized in configuration space. In terms of those functions, a localized description of polarization can be presented 
for the correlated case. Moreover, they establish the link between the Berry phase theory of polarization and Kohn's 
theory of localization . 

The following comments should be made about the above equations: unlike Eq. H, which involves an integral 
over A, Eqs. |ll] and ^2] only depend on the endpoints, A = and A = 1. The tradeoff is that whereas the former 
gives the exact change in polarization along the path, the latter give it only modulo a "quantum" [p^|-p^ (this 
"quantum of polarization" is discussed in Appendix O). A related aspect is the behavior of the equations under gauge 
transformations: unlike Eqs. |^ and |^, Eqs. |0 and |l^ are not completely gauge-invariant. As mentioned previously, 
they were obtained by assuming dxQ(k, \;G) — 0, and therefore the resulting APd is only invariant (modulo the 
quantum) under transformations which preserve the condition 0(k, A = 1; G) = 0(k, A = 0; G). Moreover, in order 
to be able to interpret Eq. |l^ at a single A as the electronic polarization |^,^, one has to impose the stronger 
condition 6(k, A; G) = 0, or — J^^'- Gauges that obey this condition are known as "periodic gauges" p2| . 

It has been pointed out by Aligia that this analysis needs to be modified in cases where there is a fractional 
number of electrons per primitive cell. The idea of a "periodic gau ge" needs to be extended to relate wavefunctions 
separated by multiples of the smallest G, and the integral in Eq. 12 needs to be extended accordingly. For the sake of 
simplicity, in the main text we will assume integer filling, and in Appendix ^ we indicate how to modify the formulas 
in order to deal with fractional filling. 

We note that in its present form the Berry phase theo ry of pola rizatio n only gives the average of the quantum 
distribution of the macroscopic polarization. In Sections IIIC and [II D we will extended it to deal with the full 



distribution. In particular, its quadratic spread will turn out to be a very relevant qua ntity, since it is intimately 
related to electronic localization. To our knowledge this was first discussed by Kudinov |2^, who proposed to measure 
the degree of localization in insulators precisely in terms of the mean-square quantum fluctuation of the ground state 
polarization. Kudinov proposed a criterion to differentiate between insulators and conductors based on the scaling 

with sample volume V of the mean square quantum fluctuation of the net dipole moment, ^Ad ^ = (^(P 

here d is the dipole moment operator, (...) means the expectation value over the ground state, and Ad = d — (d 
Using the fluctuation-dissipation theorem ||30|-|32|, Kudinov related this quantity to the optical conductivity: 



^—-L = - / du;lmxu{cj), (13) 

y Jo 

where x('^) is the electric susceptibility tensor and i = x,y,z. Using Eq. nSl Kudinov showed that as V ^ oo, 

Ad \ /V remains finite in insulators, whereas it diverges in conductors. Only finite systems with open boundary 

conditions were considered, and thus the issue of how to deal with the polarization in extended systems was not 
addressed. Our formulation leads to similar expressions, but with carefully defined bulk quantities which have a 
well-defined thermodynamic limit. 

The ideas from the Berry phase theory of polarization have recently been extended in order to address the problem 
of localization. This effort was initiated in Ref. p3| , where non- interacting electron systems with a band gap were 
considered. For such systems it is natural to attempt to quantify the degree of localization of the electrons in terms 
of the spread of the occupied Wannier functions. Marzari and Vanderbilt proposed to measure that spread via the 
quantity 



M 



[{r'),.-{r)l], (14) 



where (•.•)„ means the expectation value over the n"' occupied Wannier function in the unit cell (whose total number 
M equals the number of filled bands.) Since the electronic polarization is given by the sum of the centers of charge 
of the occupied Wannier functions [|l2|,^ , this expression is very appealing in its interpretation as the spread of the 
charge distribution of the Wannier functions. It should be noted however, that unlike the sum of the centers of charge, 
the sum of the quadratic spreads is not invariant under gauge transformations of the Wannier functions , and so 
f2 cannot be used directly as a measure of any physical quantity. Nevertheless, Marzari and Vanderbilt were able to 
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decompose it into a sum of two positive terms: a gauge-invariant part, Qi, plus a gauge-dependent term, which they 
minimized to obtain maximally-localized Wannier functions. 

We will show that fii gives the mean-square fluctuation of the bulk polarization, thus obeying a relation analogous 
to Eq. ^ (see Appendix ^). Moreover, in the same way that Q measures the spread of the Wannier functions, 
f2i measures the spread of Kohn's functions ^'mISPI, which can be interpreted as maximally-localized "many-body 
Wannier functions" ; this is discussed in Section |vrB[ and Appendix 

In Ref. it was shown that fli be rewritten as 

= T^TTs / ^kTrg(k), (15) 
[2n) Jbz 

where Tr denotes the trace, v is the volume of the unit cell, the integral is over the Brillouin zone, and g{k.) is the 
tensor 



M 



M M 



(16) 



where u„k is the cell-periodic part of the Bloch function. This tenso r is a metric which can be used to determine the 
"quantum distance" along a given path in in k-space [^3| . In Section III D we generalize this tensor to the many-body 
case, and in Section |^we relate it to the measurable polarization fluctuations via the fluctuation-dissipation relation. 

All of the above expressions for the polarization and localization involve integrals over k. More recently, alternative 
expressions have been proposed which use only periodic boundary conditions (k = 0) ||l^ , [l7| , ^^ . These are sometimes 
called "single-point" formulas. The basic quantity in this formulation is, in one dimension (ID), 



*k=o) , (17) 

where 5'k=o is the ground state many-body wavefunction obeying periodic boundary conditions over a cell of length 
L with N electrons, and as before X = X^il^i ^i- Resta [0 showed that in the thermodynamic limit the electronic 
polarization is given by 



ZN 



k=0 



^i{2w/L)X 



Pel = lim — Imlnziv. 



(18) 



However, the nature of the approximations involved at finite L, and the precise relation between Eqs. 12 and as a 
function of the size of the system, were not clarified; this is a matter of crucial importance for the usefulness of the 



expressions in practical calculations, and is discussed in Section |VII C 
Resta and Sorella 



proposed to measure the localization length in ID insulators as [p5[ 



^ = VP/(27rno), (19) 

where uq = N/ L and 

V = - lim iVlnlzArl^ (20) 

They showed that if the electrons are uncorrelated, then for insulators ^ is simply related to the ID version of ili 
defined in Eq. ^ (^^ = fli/M), and therefore ^ is finite, whereas for metals it diverges (even before taking the limit 
N —^ oo). Then they proposed that a similar behavior should occur when the electrons are correlated (with the 
difference that in general ^ diverges for correlated conductors only after taking t he lim it). 



In the present work we generalize Eq. |19| to many dimensions (see Section VII B), and give an explicit many- 
body derivation that, similarly to fli in the uncorrelated case, measures the polarization fluctuations in correlated 
systems. Similarly to the continuum formulas involving an average over twisted boundary conditions |28[ , these 
formulas require modification when there is a non-integer number of electrons per cell [p6[ (see Appendix |D|). 
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III. GENERATING FUNCTION FORMALISM 



A. Definitions 

Generating functions play a central role in the theory of statistics |^ , and have been applied to many problems in 
physics |3^j39[| . Loosely speaking, a generating function of a distribution is some function which yields, upon successive 
differentiation, the moments of the distribution, or some combination thereof. Two kinds of generating functions will 
be of interest to us: the characteristic function Cx(q:), and its logarithm, the cumulant generating function. If X is 
a vector of d variables Xi, with a normalized joint probability distribution function p{Xi, ...,Xii) — the 

characteristic function is defined as 



Cx(a)= / e-'"'^p{X,,...,Xd)dX,...dXa = {e-^'^-^) , (21) 

J — oo 

where a • X = ^t^i pQ| - The d-dimensional moments can be extracted directly from Cx(ck): 

{X-\..X--)^^-^2^, „".Cx(a) , (22) 

1 d CX — 

where n — X^iLi ^i- The cumulants are obtained in a similar way from lnCx(c*:): 

= ^.MlnCx(a) , (23) 

1 ■■■ d CX — 

where, following the notation of Ref. |Q, {■■■)^ denotes the cumulant average, which in general is different from the 
simple average (...) associated with the moments: 

(Xf)^ ^ (Xf) {X,f 
{X^X,)^ = {X,X,) - {X,) (X,) . (24) 

An important property of cumulants is that they can be explicitly represented solely in terms of the lower moments, 
and vice-versa. More precisely - and this is very relevant for what follows - for n > 1 they can be expressed in terms 
of the central moments (AX™^ AX™^ AX™^), where A.Xi = Xi ~ (Xi) and mi + m2 + ma < n, and thus they are 
independent of the mean (Xi). Moreover, provided that the characteristic function exists, the set of all the moments 
or cumulants completely determines the distribution. 



B. Polarization distribution in finite systems 



Let us consider a neutral 3D system of finite volume V containing N electrons and A^n nuclei. The dipole moment 



operator is d = qcH- + Qn^n, where X = ^t=i^i and X,i = X]i=i^i'- The average dipole moment of the system 
^I^ d ^l') = J rp{r)dr. Since the center of mass X/A^ of the electrons is not perfectly localized, ^' is not an 



eigenstate of X. Therefore the cartesian components of the dipole moment undergo quantum fluctuations, having 
a joint probability distribution p{d) dictated by the ground state wavefunction. For simplicity we will assume that 
the nuclei can be treated classically as "clamped" point charges; then they only contribute to the average of the 
distribution, and the quantum fluctuations come solely from the electrons. Hence, in what follows we will neglect 
the nuclear contribution, focusing on the distribution of the electronic center of mass. If 5* is the many-electron 
wavefunction (parametrized by the nuclear coordinates) with normalization = 1, that distribution is given by 



p(X) = (5(X-X) 
Similarly, for a given component, say, Xi, the distribution is 



(25) 
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p{Xi) = / piXi,X2,X3)dX2dX3 = SiXi - Xi) 

J — oo 

The characteristic function Cx(ck) is obtained by substituting Eq. into Eq. 



5- 



(26) 



Cx(a) - (* 



(27) 



It is clear that if we define [X]^^ X'^'^ X^'"' 



X^'X^^X!^'' 



we find, using Eq. H, (X{'iX2""X^'-^ 



X"^X2'^X3^) The electronic polarization operator is [Pc\] — QeXi/V, and the moments of its distribution 

are given by 



where n = ni + n2 + n^. 



((Pei)-(Fei)r(^ci)r> - (I)" {xrx2"'x 



* "3 

3 



(28) 



C. Polarization distribution in extended systems 

In the case of a finite system the characteristic function Cx(q^) was introduced as a purely formal device for 
obtaining the moments of the distribution, since in practice it was completely equivalent to a direct evaluation of the 
moments. In the case of an extended system the situation is rather different since, as discussed in the Introduction, the 
very definition of polarization needs to be reexamined, and a naive generalization of the direct method of calculating 
the moments does not apply. 

Again let us consider a system with N electrons in a volume V. As in previous work described in Section |l[ it 
will be convenient to use twisted boundary conditions (Eq. 0). The difficulties in defining the macroscopic bulk 

- n 1-1 

polarization can be seen from the fact that Xi is not a valid operator in the Hilbert space defined by Eq. ^: indeed, 
if l^'k > is a vector in that space, Xi j^fk > is not (it is not even normalizable), and thus (Xi ) ~ {'i'k Xi 



is ill-defined. We will now show that although the operator X is ill-defined, one can nevertheless define a meaningful 
joint probability distribution p(X) for the electronic center of mass in an extended insulating system. This will have 
the same physical interpretation as Eq. |2^ for a finite system, i.e., the moments of the variables {Pci)i = QeXi/V are 
the moments of the distribution of the electronic polarization: 



((Pel)r(Pcl)2^(Pel)r) = (I)" (^rXrXj^) . (29) 

Our main interest will be in the average (Xi) and in the quadratic spread (AX^^) which, according to Ref. [ p9| , 
measures the electronic localization. More precisely, we will define the localization length ^i along the i-th direction 
as 



lim ^(AXf). (30) 

The desired distribution p(X) can be obtained via a proper generalization of the characteristic function Cx(ck), 
which for a finite system was given by Eq. We start by introducing the quantity C(k, a), which is a precursor to 
the characteristic function for the extended system: 



C(k,a) 



,-iaX 



(31) 



where in the last equality we used Eq. ^. In the first form the similarity to Eq. ^ is evident p| ; notice, however, 
that this is not an expectation value over a single ground state \l/k, since the boundary conditions on the bra (k) 
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and on the ket (k + a) are different. That is required in order to compensate for the shift by —a in the boundary 
condition on the ket caused by the operator e~*° -^; in this way the states |^k) and 'l'k(a) ) = e"*" -'^ |^k+a) obey 



the same boundary conditions even for at ^ 2Trn/Li, so that their dotproduct ^^'k ^k(cK) y does not vanish and 

C(k, a) can be chosen to be a differentiable function (we are assuming that the ground state insulating wavefunction 
^'k is non-degenerate and is separated in energy from the excited states by a finite gap). 

That C(k, a.) is not yet the characteristic function for the extended system can be seen from the fact that the first 
"moment" that it generates is not gauge-invariant, as required from any physical quantity If C(k, a) were the 
characteristic function, the first moment would be, according to Eq. {Xi)^^ = ''■daiCCk, a)\^^Q = «($k|9fc.$k)- 
But this is nothing other than the Berry connection, which is a gauge-dependent quantity (see, for instance, Ref. 



D. Gauge-invariant cumulants 

Gauge-invariant moments and cumulants for the extended system can be obtained from the following cumulant 
generating function, which is the central quantity in this work: 



In C(a) = — ^ [ dklnC(k,a), 
(27r) J 



(32) 



where the average is over all twisted boundary conditions (for special cases where the number of electrons per primitive 
cell is not an integer, see Appendix ^). The cumulants are then obtained in the same way as in Eq. 



{X-\..X-^)^^^-^2.,.AnCia] 



(33) 



and likewise for the moments. Before we continue, we point out that C{a) is in general different from the function 



V 

Cw(q) = ^3 / rfkC(k, a). 



(34) 



This is also a characteristic function, but of a d ifferent, gauge-dependent distribution, whose interpretation and 
relation to C(a) will be discussed in Section VI B| . 

As mentioned in Section |l| Eq. |l2| is most easily interpreted when using a periodic gauge: ^k+G = ^k- Since 
we will recover that equation starting from In C(q:), a periodic gauge will be assumed in what follows; thus what is 
meant here by gauge-invariance is invariance under transformations which preserve the phase in Eq. The most 
general (/^(k) in Eq. ^ which complies with this requirement is [T^Ji^ : 



ip{k) = /3(k) - k • R, 



(35) 



where R is an arbitrary lattice vector and /3(k -I- G) — /3(k). The gauge-invariance of the cumulants generated by 
In C(q:) can now be seen as follows: according to Eqs. |9[ |2^, and |2|, the cumulants change under a general gauge 
transformation as 



{X^'X^^X. 



«3\ 

3 / 



{X^^XJl^XS-)^ + t 



n+1 



V 



dkd? 



.fc"3^(k). 



(36) 



Substitution of Eq. 35 into Eq. |3^ shows that for n = 1 the cumulants change as (X) (X) + R, whereas for n > 1 
they remain unchanged. The change in (X) but not in the higher cumulants (which do not depend on the mean) 
indicates a rigid shift by R of the whole distribution p(X.), and is related to the quantimr of polarization [|l2|,0 (see 
Appendix 



The cumulants are also real, as expected: defining /„i„2„3 



(k) =i«5«„, 



a-, Qo 



ln($k|*k+a; 



q;=0 



and /3 



-Q!, we 



obtain /„i„2«3(-k) = (-l)"i"(9^„i_^r.2^r.3 In (^-kl'S-Ck+zs)) L^g- Choosing a gauge such that $_k = *k' which can 
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always be done due to time-reversal symmetry Q, leads to ($_k |*l'-(k+/3) ) = {^k+fsl ^k), and thus /„i„2„3(— k) = 
[/nmanaCk)]*- That Completes the proof, since according to Eqs. and (^"^Xj^Xg^)^ is given by the 

average over k of /„j^„2„3(k). (Since gauge-invariance was previously established, the choice $_k = 'i'k can be made 
without loss of generality.) 

Let us now look at explicit expressions for the first few cumulants. Combining Eqs. ^ and ^and taking the first 
derivative, we find 



V 

lnC(a)|„^o = —3 
{2n) 



dk($k|9fc.*k). 



(37) 



Together with Eq. this gives 



{X,)^ = -^ j dk{^^\dM, (38) 

which, with the help of Eq. |2^, is seen to be precisely the Berry phase expression for the average electronic polarization 
(Eq. p^ ). Our formulation in terms of C(a) therefore agrees with the Berry phase theory of polarization |T2|-p^. It 
is however more general, since it also provides the higher moments: for n = 2, similar steps as before lead to 



V 



$1. 



5fe.fc,*k)-(<i>k|9fe,$k)(<i>k|5,,.<i>k) 



(39) 



Integrating the first term by parts (using the periodic gauge condition) and noting that ($k |9fc.$k) = — (c'fc^'I'kl "^k), 
this becomes {X,Xj)^ = V/{8tt^) J dkT,j{k), where 



T,,(k) = (%$k|9fc,$k) - $k) ($k |5fc^.$k) (40) 

is the gauge-invariant quantum geometric tensor [see also Eq. C9 of Ref. |^^]. The real part of Tij(k) is the 
metric tensor Gij (k) first introduced by Provost and Vallee : 



G,,{k) = Re(9fc,$k l^fc^.^k) - (%$k| $k) (*k |afc^.$k) , 



(41) 



where the second term is automatically real. Gy (k) is the many-body analogue of the tensor gij (k) defined in Eq. 
|l6|. Using Eq. 24 and the fact that the cumulants are real, (XiXj) can be rewritten as 



{X,X,) - (X,) {X,) = J dkG.,(k), (42) 

which becomes, after taking the trace, the many-body counterpart of Eq. |l^. The above equation establishes 
the physical interpretation of the k-averaged metric tensor as the mean-square fluctuation of the macroscopic bulk 
polarization (see also the next Section). The connection between such "quantum metrics" and quantum fluctuations 
was pointed out already in Ref. ||4^ . 



IV. FLUCTUATION-DISSIPATION RELATION 



We have seen that the first cumulant agrees with the Berry phase expression for the average macroscopic polarization. 
Here we will show that the expression for the second cumulant, Eq. is consistent with the fluctuation-dissipation 
relation pC|-p^ between the bulk polarization fluctuations and the optical conductivity lmxij{tj) = (l/aj)Recrt., (tj). 

For a finite system with open boundary conditions, which is the case discussed by Kudinov p9[ , the Kubo- 
Greenwood formula for the conductivity can be written in terms of the off-diagonal position matrix elements. 
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Rea"(a;) 



mlfiLoV 



(43) 



where huj„in = ^'m ^ Alternatively, the well-known relation 



-P«m = «meW„„X;„ (to ^ n) 



(44) 



can be used to rewrite Eq. ^ in terms of momentum matrix elements. Whereas for a finite system the two formulas 
are interchangeable, for an extended system the position matrix elements become ill-defined, and therefore only the 
latter form remains valid. 

At this point it is convenient to introduce the notation 



(45) 



The single-body analogue of such quantities is discussed in Refs. p3^,E6[. Our motivation for introducing them is 



the following: if P = P(*) ^^"^ many-body momentum operator for the extended system, and Pnm,k. 



^jnk), in the case of a non-relativistic hamiltonian without spin-orbit coupling Eqs. ^ and ^ lead to 



pi 



ZTOea;„m(k)X 



nm,k 



(to 7^ n), 



(46) 



which is formally identical to Eq. Notice, however, that we are now dealing with an extended system with twisted 



boundary conditions, for which the proper position matrix elements, (5'„k 
relation shows that they should be replaced by the quantities X, 



mk , 



are ill-defined; the above 



Tim,k' 



Substituting Eq. Wfl into the Kubo-Greenwood 
formula in the form valid for an extended system, we are left with an expression for Keafj^{uj) formally identical to 
Eq. H, with replaced by and w„„ replaced by w„„(k). 

Let us now specialize to the ground state {n = 0), assuming that it is non-degenerate; then Lu„ia(}i) > 0, and we 
obtain 



r t^Re<(..) = ^Re {e„>o ^o™,k^io,k} 

■ {(9A:,<i>0k I (E™>ol ^rr^i.) ($™k 1)1 4,$0k)} 



(47) 



Using the completeness relation to eliminate the excited states on the r.h.s, and comparing with Eq. we find 
Ti,{k) = Em>o^UAM- Together with G,, (k) = Rer,, (k), this yields 

^°°^Re<(c.) = ^G.,(k). (48) 
Averaging over all twisted boundary conditions and using Eqs. p3 and B2, we arrive at the desired relation: 



[ — Recrij(w), 
Jo ^ 



(49) 



where no = N/V and crij{uj) — (F/Stt^) J d\i.a^^{uj). Eq. ^ is precisely the fluctuation-dissipation relation for more 
than one variable [ ^J32| at T = (compare with Eq. [l^ ). 



V. LOCALIZATION LENGTH 
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A. Relation to the conductivity 

Here we will generalize to extended systems Kudinov's analysis of the dipole moment fluctuations as a way to 
distinghish insulators from conductors [E9j , and discuss it in terms of a localization length. It is convenient for that 
purpose to classify solids into three categories, according to the low-frequency behavior of the conductivity at T = 
as V ^ oo: 



lim^j^o Re(T(a;) = insulators 

Rea{uj) = (27rgg/?i^)£>(5(a>) +Ref7rcg(w) ideal conductors (50) 
lim^^^o Re(T(w) = ctq non-ideal conductors 

Insulators are characterized by a vanishing dc conductivity, in contrast to conductors. The singular contribution 
occurs in ideal conductors - those without any scattering mechanism - and D is called the "Drude 
weight" or "charge stiffness" p9| , ^ . If there is scattering, the S function peak is smeared out to a Lorentzian, so 
that D = and the dc conductivity of non- ideal conductors, such as disordered metallic alloys, has a finite value (Tq! 
(Trog(w) is the regular finite-frequency part of a{u)) in perfect conductors. 

The only possible divergence of the integral on the r.h.s. of Eq. ^ is around cj = H, because as ^ oo 
luiXiiito) ''^ (^~^ Substituting Eq. ^ into Eq. ^ we see that the quantity in Eq. 30[is finite for insulators 
and diverges for conductors. By the same token, using the zero- frequency limit of the Kramers-Kronig relation 
one finds the familiar result that the static susceptibility is finite for insulators and divergent for conductors. Hence, 
assuming Eq. |5^, at T = the following three conditions are equivalent: (i) Rexfj(O) is finite, (ii) lim^j^o Ke(Tii{uj) — 0, 
and (iii) is finite. 

The quantity has the dimensions of length. Since it is finite for insulators and infinite for conductors, it is natural 
to interpret it as an electronic localization length along the i-th direction. According to Eq. ^ for extended systems 
it can be written as = lim^v^oo Ci(-^)j where ^i{N) is given in Table I in terms of the cumulant generating function. 
This formula has a striking similarity to the Drude weight formula derived by Kohn in terms of the total energy E(\i) 
[p|, ^7|j4S| ] , given also in Table I: both are second derivatives of some quantity with respect to the twisted boundary 
conditions (in the case of (N) the twisting of the boundary condition k is followed by an averaging over all k, hence 
the parameter a instead of k). 

It is clear from Table I that unlike the localization length, the Drude weight does not provide a universal criterion 
to discern insulators from conductors However, the combination of the two quantities in principle enables us to 

distinguish between the three categories. In the same way that the Drude weight measures the "degree of conductivity" 
of an ideal conductor, 1/^^ measures the degree of localization of the electrons in an insulator. Insofar as localization 
- in a properly defined sense - is an essential property of the insulating state [^|j|] , this can be viewed as a meaningful 
measure of the "degree of insulation", one which is expected to apply to all types of insulators. 

Although we have managed to express in terms of the measurable optical conductivity via the fluctuation- 
dissipation relation, it is not yet clear how it relates to the notion of localization put forward by Kohn, in terms of the 



localization properties of the insulating wavefunction in configuration space. That will be discussed in Section VI B 



B. Relation to the optical gap 

In an insulator optical absorption starts at a threshold energy i?g, below which Keaii{uj) — (we are neglecting 
phonon- assisted transitions, so that the gap Eg is the minimum gap for optical transitions). From this it follows that 



Jo ^ Eg Jq 



H- (51) 



With the help of the sum rule for oscillator strengths ||l|,|^,^ dLjReaii{uj) — (l/8)Wp (wp is the plasma frequency), 
together with Eqs. |3^ and we conclude that 

znipEg 

This inequality shows that the polarization fluctuations are controlled by the optical gap, lending support to the 
intuitive notion that the larger the gap, the more localized the electrons are. It strongly resembles an inequality 
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previously derived by Kivelson [ |5l[ | for non-interacting electrons in ID, where is replaced by the quadratic spread 
of properly chosen Wannier functions. As discussed in Appendix |a[ in fact in ID equals the average spread of the 
maximally-localized Wannier functions. 



VI. "MANY-BODY WANNIER FUNCTIONS 



The expression for the cumulant generating function, Eq. |3^, involves an average over all twisted boundary condi- 
tions, which was introduced in a somewhat ad hoc manner in order to render the resulting distribution gauge-invariant 
(modulo a rigid shift by a quantum) . We will now shed some light on the physical significance of the averaging pro- 
cedure, by showing how it can be rationalized in terms of the notion of localization in insulators developed by Kohn 
This will be achieved by introducing "many-body Wannier functions", and will allow us to tie together the 
Berry phase theory of polarization and Kohn's theory of localization. 



A. Finite system: Localization in real space 

It is instructive to start by discussing the case of a finite system of linear dimensions ~ a (e.g., a molecule with 
N electrons). Instead of the usual open boundary conditions, we can choose to impose periodic boundary conditions 
on its wavefunction (k = in Eq. ^), choosing L » a. The resulting A^-electron wavefunction is periodic in 
configuration space, as depicted schematically in Fig. 1 for ID and A'^ = 2; it is a sum of partial wavefunctions which 
only differ from one another by a translation of the coordinates. 



*k=o(xi, ...,Xjv) =^...^^'mi...m„(xi,...,XAr), (53) 
mi m^v 

where the {nii} label the partial wavefunctions. These are localized in geometrically equivalent regions 7?.mi...mjv 
configuration space (shaded regions in Fig. 1), which for L » a are essentially non-overlapping: 



«'mi...m„^'rn;...m;^ = 0' for (nil , . . . , HiAT ) 7^ (mi , . . . , m^y) , (54) 

where, using the notation of Refs. the symbol = denotes equality apart from exponential small corrections which 
vanish in a manner such as e~^/^, where ^ is a localization length (in this example ^ a). 

Next we switch from periodic to twisted boundary conditions (k ^ in Eq. From the localization properties it 
follows [see Section 2 of Ref. and Section VI of Ref. that the periodic part of the wavefunction, $k = e~^^ '^'^\i, 
can be written with an exponentially small error as 



$k(xi,...,xjv) =e *''<^$k=o(xi,...,xjv), (55) 
where we have introduced the quantity 

Q — / Ei=l(Xi ~ in 7?.mi...mjv fr^Q'j 

\ F(xi, ...,XAr) outside aU7?.mi...m„ ' 

where = Lrxii, and F is a largely arbitrary, periodic function, which joins smoothly with the values of Q at the 
boundaries of the regions 7?.nii...mjv Let us now look at the precursor characteristic function C(k, a) for such a finite 
system obeying twisted boundary conditions: substituting the previous expressions into Eq. 0, we find 



C(k,a) = ($k=o 



$k=o ) = ( * 



^o...o 



0...0 



(57) 



which has the form of Eq. obtained using open boundary conditions. It is clear that as a result of localization of 
the electrons in real space the choice of boundary conditions is immaterial. In particular, C(k, a) becomes independent 
of k, and thus from Eqs. 32 and we conclude that 
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C{k,a) = Cw{a) = C{a). 



(58) 



Hence, averaging over boundary conditions becomes redundant, and C(k, a) is already the correct, gauge-invariant 
characteristic function. Finally, it is interesting to note that Eqs. ^ and |5^ can be viewed as a particular realization 
of Eqs. 2.22 and 2.23 of Ref. which were used to relate the quantum metric to the quantum fluctuations. 

B. Insulators: Localization in configuration space 

Let us now consider an extended system, such as a crystal with twisted boundary conditions over a large volume 
V. Although in general the charge density will be delocalized in real space, Kohn has argued that if the system is 
insulating the wavefunction is localized in configuration space . Kohn's notion of localization is weaker than the 
one implied by Eq. ^ and Fig. 1; however, as we will see, it is sufficient to recover Eq. |5^, after making a judicious 
choice of gauge. 

We start by introducing a localized description of the many-electron insulating wavefunction in configuration space. 
In the non-interacting case, a localized one-electron description (in real space) is provided by the Wannier functions 
|p3| , p6|j43| . The Wannier function |Rn) associated with band n and centered around the unit cell labelled by the 
lattice vector R is related to the Bloch functions by the following unitary transformation: 

|Rn) = ^^e-'-^|V'„k), (59) 

where a periodic gauge is assumed, and the sum is over a uniform grid of Nc points in the Brillouin zone p2[ . Due 
to the discretization of the integral over k, |Rn) is actually periodic in a large cell of volume N^v (v is the volume 
of the unit cell). The many-body analogue of a periodic Wannier function ("many-body Wannier function") can be 
defined in a similar way in terms of the many-body wavefunction (again assuming a periodic gauge) : 

\W^) ^ ^J2e-^^-^^ \^^) , (60) 

SO that, using Eq. ^, 

|$k) --i^y^e-*''-(^-^") |M^m). (61) 



Here Rm = where the Mi are integers, and the vectors {Li,L2,L3} define the volume V containing 

N electrons (notice that V is the volume of a supercell, typically a large multiple of the unit cell volume v). $k is 
periodic and 5'k obeys k boundary conditions over the volume in configuration space, and both ^'k and Wm 
are periodic over the volume (NcV)^ . The normalization conventions are the following: $k is normalized to one 
over a volume , and Wm is normalized to one over a volume (NcV)^ . Similarly to the |R?i), the |H/m) form an 
orthonormal set. 

From the same type of general considerations which are used to show that Wannier functions can be chosen localized 
p6[ , it follows that we can choose Wm localized in the variable X = X^i^i with a distribution 



Pw(X) = (Wm S{X - X) Wm) ■ (62) 

Substituting Eq. ^ in Eq. averaging over k, and comparing with Eq. |3j, after discretizing the integral over k, 
yields: 



Cw(a) ^^Y.^{-k,a)^(wM e-'°-^ Wm) , (63) 



Nc 
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which shows that the function C^^{a.) introduced in Section HID is the characteristic function of the distribution 
(X) . It is straightforward to check that the first moment of this distribution equals the first moment of the gauge- 
invariant polarization distribution p(X) generated by C{a); this is similar to the independent-electron case, where 
the average polarization is given by the sum of the centers of charge of the gauge-dependent Wannier functions . 
Likewise, the "many-body Wannier functions" W-m are gauge-dependent, and as a result so are the higher cumulants 
of pw(X), in particular the quadratric spread. Notice that the gauge-dependence of the Wm. implies that in general 
they overlap with one another , and only in certain gauges are they localized to the point of being essentially non- 
overlapping (provided that the system is insulating). This is different from the usual single-body Wannier functions, 
which in general remain overlapping in any gauge (even though they are orthogonal); in the many-body case we are 
free to choose the volume V large enough so that with a judicious choice of gauge the Wm become non-overlapping 
in the high-dimensional configuration space. 

We are now in a position to provide the link between the present formalism (and, therefore, the Berry phase theory 
of polarization) and Kohn's theory of localization ||,^. A major result of Kohn's work is the conjecture that a 
general many-body insulating wavefunction ^k=o breaks up into a sum of non-overlapping parts localized in 
disconnected regions TZm in configuration space (Eq. 0). Upon inspection, one finds that the pre sent "many-body 
Wannier functions" are nothing other than Kohn's functions, Wm. = v'^^'m [compare Eq. pl^ with Eq. 6.1 of 
Ref. 1^], except that Kohn only considered gauges where they are localized in such a way that the overlap becomes 
exponentially small: 



*M(xi,...,XAr)«'M'(xi,...,xjv) = forM^M, (64) 

so that the ^'m are uniquely defined apart from exponentially small variations and a global phase. Transposing the 
language of Ref. to a many-body framework, Kohn's functions ^'m can be viewed as the maximally-localized 
"many-body Wannier functions" Wm (see Appendix |b]) . 
In such "non-overlapping gauges", the precusor characteristic function, Eq. becomes, using Eqs. |6l| and |6^, 



'M 



^'m) ("non-overlapping gauge"), 



(65) 



which is independent of k, since the k-dependence of C(k, a.) in a general gauge arises from the cross terms between 
different Wm- Therefore in such gauges we recover Eq. |58|, which for a finite system obeying twisted boundary 
conditions was valid in any gauge: 



C(k, a) = Cw(a) = C{a) ("non-overlapping gauge"). 
Finally, from the gauge-invariance of C(q:) we conclude that in any gauge 



(66) 



'M 



M 



(67) 



which is the desired relation. It should be stressed that in deriving this equation we assumed Kohn's conjecture 
regarding the existence of non-overlapping many-body Wannier functions in the insulating state. 

We interpret Eq. ^ as follows: in each region TZm in the configuration space of the variables {xi} (Fig. 2), the 



variable X = X]t=i ^* takes on a range of values with a distribution generated by C(a), 



p(X) =iV,(*M <5(X-X) ^ 



'M 



(68) 



which, together with Eq. |2^, gives the bulk polarization distribution. The previous equation can therefore be viewed 
as the generalization for extended insulating systems of Eq. |2^, which gives the distribution of the electronic center 
of mass for finite systems. 

In Fig. 3 is represented the distribution p{Xi) along the i-th direction. According to Eq. |3^, the width of this 
distribution is ^/N^i(N). The solid lines describe what happens in a ID insulator, for which the width of each of the 
peaks p{X) labelled by M is \/N^{N) oc and the distance between the centers of consecutive peaks is L, so that 
for large L they are well-separated. In the case of a 3D insulator where all the linear dimensions are similar {V ^ L^), 
the peaks p{Xi) overlap for large L, since their width is ^/N^i{N) oc L^^^ (dashed lines in Fig. 3.) The important 
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observation is that, according to Eqs. ^ and even if they are overlapping, the distributions p(X) are well-defined 
"projections" into real space of essentially disconnected distributions I^'mI^ in configuration space js^ . 

At this point we shall reconsider the problem of how to define a meaningful many-body position operator for the 
electrons in extended systems, an issue which has been recently discussed in Ref. We recall that the source 

of the problem is that in such systems the usual center-of-mass operator X is ill-defined (see Section III C ) . Kohn 
proposed the following operator as the substitute for X for taking expectation values over the ground state 
wavefunction of an extended insulator. 



Q 



X — R.M in T^M 

F(xi, ...,X7v) outside all TZm ' 



(69) 



where F is defined in a similar way as in Eq. |5^. It is straightforward to check that the moments calculated from Q 
coincide with those derived from the characteristic function: with the help of Eqs. 0, and |6^, we find 



'M=0 



*k=0 



o) =^"a;VC(a) 



a=0 



(70) 



The connection with Ref. [O] will be made in Section VII B 



C. Insensitivity of bulk properties to the boundary conditions 

We now have an adequate framework for discussing the implications of localization for the dependence on boundary 
conditions of the bulk properties of insulators. As a first example, consider the adiabatic current Since it is 

gauge-invariant, we may evaluate it in any gauge, in particular in a "non-overlapping gauge" . Substituting Eq. |6l| 
into Eq. I?] we find, because of the exponential decrease of the overlap as the size L is increased: 



M 



X 



(71) 



where for the last equality Eqs. ^ and ^ were used. Thus for a large system size j'^'' is essentially k-independent 
(this was demonstrated in Ref. using a different reasoning, for insulators with an energy gap). 

As a second example, let us look at the quantum geometric tensor TijCk) and its real part, the metric tensor (k). 
Evaluating Eq. ^ in a "non-overlapping gauge" with the help of Eq. we find Tij(k) = (XiXj) — (Xi) (Xj), which 
again is essentially independent of k for large systems (it is also real, so that T!y (k) = (k)). Thus we conclude that 
the k-independent metric tensor gives the mean-square fluctuation of the polarization. This is a stronger statement 
than the one made at the end of Section HID , which pertained to the k-averaged metric (by contrast, the single- 
electron metric tensor gij(k) given by Eq. |16| is in general k-dependent, since the one-electron Wannier functions 
remain overlapping even when they are maximally-localized; as a consequence, the quadratic polarization fluctuations 
are related to its average over the Brillouin zone, as shown in Appendix ^.) 

These two examples illustrate an important point: the physically measurable ground state quantities, such as the 
polarization current and the polarization fluctuations, are well-defined for any single boundary condition k, and the 
choice of k becomes immaterial for large system sizes. It may therefore seem bizarre that in the generating function 
formalism one needs to average over all k before obtaining gauge- invariant quantities, as shown in Section HID. A 
concrete example helps to clarify this state of affairs: in the case of the second cumulant, had we taken the second 
derivative of the precursor cumulant generating function, lnC(k, a), without afterwards averaging over k, we would 
be left with a quantity (the integrand of Eq. ^), which is different from the desired one, — Ty (k), and is gauge- 



^uge-dependent quantity is gauge- invariant, and equals the k- average of 
3). A similar situation occurs with the adiabatic polarization current. 
The insensitivity of the bulk properties of insulators to the choice of boundary conditions on a large system seems 
to be a very general property. Indeed, if O is any well-defined operator acting on the ground-state of the extended 



dependent. However, the k- average of that 
—Tij(k.) = — Gy (k) (compare Eqs. pa and k 



«'k=0 



O 



system (e.g., the hamiltonian H), then from Eqs. ||, and ^ it follows that ^\E'k O ^'i 

which helps to understand why, for instance, in insulators the Drude weight (see Table I) goes to zero exponentially 
with the system size Q. This has been confirmed by numerical simulations and calculations on exactly solvable 
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models Q. More precisely, in Ref. |Q it was found that the Drude weight for the half- filled Hubbard model in ID 
scales as ~ L^/^e~^/^, which yields a localization length ^. It would be worthwhile to investigate whether ^ is the 
same as the localization length ^ used in the present work. 



VII. DISCRETIZED FORMULAS 



Until now we have dealt with a continuum formalism, in which the cumulants are obtained by differentiating 
In C(q:) at a = 0. Heuristically, this can be viewed as a measure of the k-averaged change in $k as the boundary 
condition k is twisted adiabatically (see Eqs. ^ and |3^). In numerical calculations it is often more natural to 
perform independent calculations on a discrete mesh of k-points and then use the resulting wavefunctions to estimate 
the derivatives by finite differences |]l^,^. Since the global phases of the wavefunctions on the grid are unrelated, 
useful discretized expressions must remain invariant under arbitrary changes of those phases. The derivation of such 
discretized formulas becomes particularly transparent in the present formalism, as we will now show. 



A. Formulas involving an average over k 



We will find it convenient to work with scaled coordinates: let the simulation cell be defined by the vectors 
{Li,L2,L3}, to which correspond the reciprocal lattice vectors {G^, G^, G-^} (Li-G^ — Stt^/). The scaled coordinates 

are defined by X = S'^L^, where a sum over repeated indices is implied, and similarly we have k = QiG^ and a — UG^ 
(and thus a ■ ~ 2'KtiS^). In terms of these variables, Eqs. Isil and ^ become 



C(q,t) = (*q 



-i2iTtiS^ 



*q+t) = ($q|$q+t), 



and 



lnC(t)= / dqi [ dq2 f dq3\nC{q,t) 
Jo Jo Jo 



(72) 



(73) 



so that C(t) is the joint characteristic function for the variables 2ttS'^ 
The average is 



{S') = ^ J dqdt, lnC(q,t)|,^o. 
Next we discretize the derivative in the integrand as 



(74) 



6qidt, lnC(q,t)|,^o ~ lnC(q,(5qO - lnC(q,0) ~ iImlnC(q, (5qO, (75) 

where we made use of Eq. This expression is gauge-dependent, just like its continuum counterpart, the Berry 
connection. As in the continuum case, gauge-invariance is recovered by averaging over k; that is done by choosing a 
uniform row of J+1 k-points along the direction of G' such that the endpoints are separated by G', k-y = k^ + (7/ J)G' 
with 7 = 0, J; then Sqi = 1/J, and we find: 




where is some permutation of (1,2,3), and the periodic gauge is enforced by setting <I>kj — e '"^'I'ko- The 

average polarization is then given by 



(Pel) • G' = 
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which is the many-body analogue of the discretized Berry phase formula proposed in Ref. [|2|. It is straightforward 
to verify that it is gauge-invariant. This property hinges upon (i) using a uniform row of k-points, (ii) choosing a 
periodic gauge, and (iii) having a logarithm in the expression, which in the present derivation appears quite naturally, 
coming from the cumulant generating function, Eq. 
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To compute the variance we discretize the second logarithmic derivative: 



t=0 



lnC(q, 5(li) + lnC(q, - Jq;) - 2 lnC(q, 0) = In [C(q, (5q,)C(q, -(5q,)] 



(78) 



which is gauge-dependent, similarly to Eq. |75| ||5g]. After some manipulations we obtain the following gauge-invariant, 
k-averaged formula: 



(5')- 



J 



dqidqj In 



J-i 

n (*k,i<i>k,+i) 

7=0 



The evaluation of the covariance requires taking cross derivatives, which can be discretized as 

^xyf - 2sts^ [f(^ + Sx,y + 6y) + f{x ~Sx,y~ 6y) 



f{x + Sx, y) - f{x - 5x, y) - f{x, y + 5y) - f{x, y - 5y) + 2f{x, y)] , 



leading to 



^^ntio' u':J K<fk.. l*k.+,„ )|' - inn.'ro' u''J K<fk.. I'l'k.,.^. )| 



(79) 



(80) 



where K-y = ^± + j-G^ + ^GK 

In order to convert back to cartesian coordinates we use the transformation law for a second-rank contravariant 
tensor, (XiXj)^ = Hii{S'' S"") ^Hjm, where H = {Li, L2, L3}. Eqs. |7^, |7|, and |l| are all that is needed to calculate 
the average polarization and its quadratic fluctuations, as well as the correlations along different cartesian directions 
if the cell symmetry is orthorhombic or higher, Eq. ^ is not needed for computing the fluctuations.) In Appendix 
we give the required modiflcations to deal with fractional filling. 



B. Formulas involving a single k 

We now turn to a different kind of discretized formulas, the so-called "single-point" formulas, such as Eq. In 
the case of localization they were previously given for ID only ||3^ . In the present formalism the generalization to 
higher dimensions becomes straightforward. We will arrive at the "single-point" formulas starting from the expressions 
derived in the previous Section, which involve averages over a uniform grid of k-points. The basic idea is to perform 
the many-body analogue of a "Brillouin zone folding" p6[ |. Let us start by discretizing the remaining integrals in Eqs. 
[TT] , |79| , and 81; the expression for the average polarization, for example, becomes 



V J^JJ 



J,-l ,7,-1 

-ImlnJ] n n (*k_|$k.,.,.,,) 
p=0 u=o 7=0 



(82) 



where we have used a uniform spacing along all the reciprocal lattice directions: k^^^.^ = k-|- (^/Ji)G' -I- {v/Jj)G^ + 
(7/Jz)G' and k is fixed (usually k = 0.) Next we build an ansatz wavefunction jsj] contaning N — JN electrons 
as the antisymmetrized product of the J = J1J2J3 separate A''-electron wavefunctions "fk^^.^: 



1 



:^(-l)^P*(l,...,iV), 



(83) 



17 



where the sum is over all the permutations P of the N particle coordinates, (—1)^ is the parity of the permutation, 
and ^'(1, N) = Yip^u-y ^k^„^- It is easy to verify that 5"^^ obeys k boundary conditions on the larger cell containing 
N electrons. It can be shown that for even N this leads to (compare with Eq. 16 in Ref. jsj], which deals with the 
ID case): 



k/ ' 



(84) 



where G' = / Ji is a basis vector of the reciprocal of the larger cell, and X is the center of mass position operator 
of the TV-electron system. Comparison with Eq. p3 yields 



(Pei)-G'~-4^Imln(vl/,~ 



V 



-iG'-X 



(85) 



Similarly, discretizing the integrals in Eq. ^ and using Eq. we find 

1 Ji 



{2TJ:f JtJ, 



■In 



iG -X 



(86) 



The quantity on the l.h.s. pertains to the A^-particle system, whereas those on the r.h.s. pertain to the A^-particle 
system. However, if we choose Ji = J2 = J3 and use {XiXj)^/N ~ (XiXj) /N together with the transformation 



law between cartesian and scaled coordinates, we find (^S'S^j = Ji{S'S^)^. Eq. |86|then becomes 

1 



{s-y 



■In 



-iG'-X 



where all the quantities are now explicitly written for the A^-particle system. 
The following relation can be derived along the same lines as Eq. p4: 



(87) 



g-iG^-Xg-iG'-X 



n('fk..-,i*k,„+i,.+i)=(^t 

Discretizing the remaining integral in Eq. ^ and again setting Ji — J2 — J3, we obtain, using Eqs. ^ and [88 
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(89) 



and|8^ are the desired "single-point" formulas. In Appendix ^ we modify them for the case of fractional 



Eqs. 1^ 
filling. 

In ID the above formulas become particularly simple, and we recover the results of Rcfs. and In cartesian 
coordinates {X — SL and G = 2tt/L) they are 



X 



which is equivalent to Eq. |l^, and 



— Imln /^r 
27r \ 



-i^X 



k / ' 



(90) 



■In 



^1 



(91) 



which, together with Eq. BOl gives Eq. pi 
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C. Interpretation of the "single-point" formulas 



In the previous Section we derived "single-point" formulas for the first two moments of the polarization, starting 
from discrertized expressions involving an average over boundary conditions. In order to switch between the two 
descriptions, an ansatz wavefunction given by Eq. ^Sj^was used. However, this is strictly valid only if the N 
particles are not correlated with each other. In that case becomes a Slater determinant of one-electron orbitals, 
and the procedure leading to the "single-point" formulas ("Brillouin zone folding") does not involve any further 
approximations (see Refs. and p7[|). 

For a correlated state of many particles the situation is rather different. In that context the "single-point" formulas 
were originally proposed in ID for an arbitrary correlated iV-electron wavefunction ^'^^g with periodic boundary 
conditions over a cell of size L \17 34 



such wavefunction in general will not obey Eq. |8^, of course. The derivation 
of the previous Section allows us to assess the approximations involved: these have to do with the extent to which 
a wavefunction given by Eq. ^ differs from the fully correlated wavefunction. The key quantity to consider is the 
correlation length Lcom which quantifies the range over which the particles are correlated (in addition to any long 
range order which is included in the mean field potentials seen by each particle). The fact that the correlations are 
short-ranged has been termed "nearsightedness" by Kohn in a recent paper |^7|| ; in fact, in an insulator, for example, 
one expects the longest range correlations to be of the Van der Waals type, which decay as 1/r^ in the energy and as 
in the wavefunction |39,Q. The basic assumption underlying Eq. namely that the N electrons only correlate 
in groups of at a time [34], is consistent with the principle of nearsightedness. This is the justification for applying 
the "single-point" formulas to a correlated insulating wavefunction, provided that the cell is larg e e nough. 

From this perspective it becomes clear that the two types of formulas derived in Sections VII A and VII B constitute 
different approximations to the same continuum expressions, which involve an average over all twisted boundary 
conditions. In both approaches one must choose cells with linear dimensions greater than Lcorr: in the latter approach 
one must use a single large cell with size JV in order to have the same level of accuracy as J independent calculations 
each of size V using the former approach. 

From a practical point of view, there are two possible ways to proceed: either perform several calculations with 
different twisted boundary conditions on a smaller cell and then use the formulas of Section VII A , or perform a single 

Although the two approaches are comparable 



calculation on a large cell and then use the formulas of Section VII B 



in terms of accuracy, the former approach should be more efficient, because the available many-body algorithms do 
not scale linearly with the number of particles. On the other hand, the averaging over boundary conditions may be 
more cumbersome to implement in practice. 



VIII. CONCLUSIONS 



We have presented in this work a unified theory of electronic polarization and localization in insulators. The central 
quantity in the formalism is the cumulant generating function In C(q:) defined in Eq. ^2|; it provides a systematic 
procedure for extracting from the ground state wavefunction \['k the moments of a properly defined distribution p(X) 
for the electronic center of mass in an extended insulator. In complete analogy with the case of a finite system, this 
distribution is simply related to the quantum distribution of the ground state polarization, via Eq. |29|. 

Several seemingly disparate ideas regarding electronic polarization and localization in insulators p|, ^4| , ^9[ are brought 
together quite naturally using the generating function approach. In particular, it shows the connection between the 

and Kohn's theory 
rich is defined in terms of 



Berry phase theory of polarization [ [12|jl4[ , as well as the related theory of localization ||33 
of localization in the insulating state |B|,p|. A key quantity is the localization length ^i, wt 
the experimentally measurable mean-square fiuctuation of the polarization (Eq. ^ and Table I), making contact with 
the work of Kudinov In the thermodynamic limit agrees with the localization length defined in Ref. ^J] for ID 
systems (Eq. p^ ). Furthermore, the generating function formalism also reveals a very close formal analogy between 
Kohn's localized many-electron functions and the maximally-localized one-electron Wannier functions defined 
by Marzari and V anderbilt |^ : the former can be viewed as maximally-localized "many-body Wannier functions" 
(see Section VI B and Appendix ^). Moreover, in the same way that the quadratic spread of the functions ^'m is 
a measure of the mean-square fluctuations of the bulk polarization, the gauge-invariant part of the spread of the 
Wannier functions (which in ID systems equals the spread of the maximally-localized Wannier functions) measures 
the same quantity for uncorrelated insulators. The fluctuation-dissipation relation can be used to derive an inequality 
(Eq. |5^) between the polarization fluctuations and the minimum energy gap for optical absorption in an insulator. 
The present approach also provides some extra insight into the appearence of a "quantum of polarization" in periodic 
insulators , which is seen to be related to the existence of multiple geometrically equivalent regions TZ-wi , depicted 
in Fig. 2 (see Appendix ^). 
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The localization length seems to play a role in the theory of insulators similar to that of the Drude weight in the 
theory of ideal conductors: the latter measures how "free" the "free charges" in a perfect conductor are, whereas 
the former measures how localized the "bound charges" in an insulator are. Interestingly, both quantities can be 
expressed as second derivatives with respect to th e twisted boundary conditions (see Table I), which play a crucial 
role in the formalism. As discussed in Section VI C, bulk properties of insulators arc rather insensitive to the boundary 
conditions, unlike the properties of conductors. Nevertheless, the insulating wavefunction itself is not insensitive to the 
boundary conditions [ p5| ; on the contrary, the derivatives of In C(q:), which measure the k-averaged change with k in 
$k = e^'*' "^^k as the boundary condition on is twisted, contain quantitative information about basic properties of 
the insulator: in particular, the first derivative gives the average macroscopic polarization, and the second derivative 
gives its mean-square fluctuation. Both quantities have a geometrical interpretation: the former is a Berry phase on 
a manifold of quantum states parametrized by the twisted boundary conditions (l^Jl^ (Eq. |3^), and the latter is a 
metric on the same manifold (Eq. |4^ ). 

The generating function approach also leads naturally to discretized formulas which can be used to compute the 
polarization and the localization in many-body numerical calculations in any number of dimensions. Two alternative 
kinds of expressions exist: those involving wavefunctions computed on a uniform grid of k-points (Eqs. and 
|8l| ), and those involving a single wavefunction with a fixed boundary condition k (Eqs. and^). The present 

derivation clearly shows how the two types of formulas are related to one another. 
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APPENDIX A: LOCALIZATION LENGTH AND WIDTH OF THE WANNIER FUNCTIONS 



In Section VI B it was shown that the width of the electronic center of mass distribution p{Xi), arising from the 
maximally- localized "many-body Wannier function" 'I'm, is \/N£,i. In this Appendix we work in the independent- 
electron framework, in which the usual one-electron Wannier functions are defined, and investigate the relation 
between their width and the localization length (i.e., the mean-square fluctuation of the polarization). It is clear 
that since the spread (Eq. ^4|) is gauge-dependent, it cannot relate directly to any measurable quantity. One might 
have guessed that, as happens with the "many-body Wannier functions", the physically meaningful quantity would 
be the spread of the maximally-localized Wannier functions. Building upon the results of Ref. we show here 
that the gauge-invariant part Qi of the spread of the occupied Wannier functions (Eq. [l5| ) measures the mean-square 
fluctuation of the bulk polarization. 

The proof follows from the fluctuation-dissipation theorem: for a crystalline insulator in the independent-electron 
approximation, the real part of the optical conductivity due to interband vertical transitions is given by |50[ 



M 



where / is the occupation number of states in the valence band (in spin-degenerate systems / = 2), p„m.k = 
—ih (V'nkI VlT/imk), and huj^^n — ^mk ~ ^nk IS the difference between single-particle energies. Instead of Eq. Wq we now 



have Pnm k — ^e^mn (^nkl^kUmk), and foUowing similar steps as in Section IV for the many-body case we find, using 
Eq. 0, 



-^Rea,(.) = ^y^^dk,,(k). 



(A2) 



Combining with Eqs. and |9[ using v/{fV) — AI/N, and taking the trace, we obtain the fluctuation-dissipation 
relation |59| ]: 



20 



n,^^ r ^TrRea(c.) = M lim (A3) 

Comparison with Eq. ^ shows that — X]i=i ^fi which is the multidimensional generalization of a result obtained 

in Ref. |^^. It has been found that typically Slj accounts for most of the spread ft; for instance, in semiconductors it 
usually accounts for more than 90% [Q. Thus it is justifiable to view as an estimate of the average width along the 
i-th direction of the occupied Wannier functions (in ID is actually the average spread of the maximally-localized 
ones Q). This conclusion agrees with the analysis in Appendix 1 of Ref. there, the distribution p{X) of the 
electronic center of mass of a ID non-interacting insulating system was studied, with the result that X is localized in 
an interval of width \/Wb, where "& is a length of the order of the range of the Wannier function" . It is clear from 



the results of Section VI B and the previous discussion that the gauge-invariant ^ is precisely that length b. 

It is straightforward to check that for non-interacting crystalline insulators, Eq. ^ becomes, after summing over 
all d cartesian directions, fli/AI < dh^ / {2mf,E^), where is the minimum direct gap over the Brillouin zone. Hence 
in general the inequality involves the gauge-invariant part of the spread. Since in ID this equals the average spread 
of the maximally-localized Wannier functions, we recover Kivelson's original result pO|. 



APPENDIX B: SPREAD OF THE "MANY-BODY WANNIER FUNCTIONS" 



In Section VI B we introduced the functions Wm , which play a role in the many-body theory of polarization similar 



to that of the Wannier functions in the independent-electron theory. Here we will show that this formal analogy carries 
over to considerations about the spread of those functions. As mentioned in Section ||, in the independent-electron 
framework the quadratic spread f2 of the occupied Wannier functions can be decomposed into a sum of two positive 
terms, one of them (ilj) gauge-invariant |^^; hence fij < $7 in any gauge. Following similar steps as in Ref. [ p3[ , 
here we will derive the corresponding result for the many-body functions Wm- First we define their gauge-dependent 
spread by analogy with Eq. Ill: 



X^ 



Wn 



Wo) =-diAnCw{a)\ 



(Bl) 



where C\n{ol) is given by Eq. Next we show that the role of Vli is played by (AX^) = ~^a2 ^'^^('^)|q;-o' ^^'^^S 
this derivative inside the integral in Eq. p3, discretizing that integral 152] , and comparing with Eq. B3, this becomes 



(AX2) = (Wo 



X^ 



Wo 



(B2) 



Now we use the relation 



M 



X 



Wo), 



(B3) 



which is the many-body analogue of Eq. 5 of Ref. | |33| and can be derived in the same way. Substituting into Eq. B2 
gives 



(AX^) = (Wo 



X2 


Wo)-^ (Wm 


X 


Wo) 




M 







(B4) 



and comparison with Eq. Bl yields the desired result: (AX^) < 17^^'. 

There is, however, an important difference with respect to the single-electron Wannier functions: in Ref. p3| it was 
shown that it is only in ID that the maximally-localized one-electron Wannier functions have a spread min[f2] — f^i 
whereas in higher dimensions min[S7] 
according to Kohn 



> Oi, i.e., fii is strictly a lower bound to fi. As discussed in Section VIB, 
a gauge can be chosen where the functions VFm have an exponentially small overlap with one 
another in the high-dimensional configuration space, in which case only the term M = survives in the sum on the 
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r.h.s. of Eq. B4. Therefore we conclude that for an insulator in any number of dimensions we have min[r2(^'] 
•-I 



(AX^), where for any finite size the corrections due to the exponentially small overlaps make the l.h.s. slightly larger 
than the r.h.s. Thus Kohn's functions \['m are the maximally-localized "many-body Wannier functions" H^m, since as 
y — > oo the width <^AX^^ of the gauge-invariant distribution p(X) obtained from (Eq. ^8|) becomes the minimum 
of the width of the gauge-dependent distribution pw(X) obtained from W^m (Eq. |62| ). 

APPENDIX C: THE QUANTUM OF POLARIZATION 

Here we will discuss the quantum of polarization in terms of Kohn's theory of localization |||,^ [A related 

discussion, regarding the closely related phenomenom of quantized charge transport in insulators, can be found in Ref. 
1^.] Let us consider a periodic system of volume ^ which for simplicity we take to be cubic. If the external potential 
acting on the electrons is changed adiabatically along an insulating closed path parametrized by A {H^^^^'' = H^^^°^), 
then, since the hamiltonian comes back to itself, the net effect, as far as the wavefunction is concerned, has to be 
either (i) each returns to itself (modulo a global phase, the same for all M), or (ii) apart from a global phase, 
there is a rigid translation of all the ^'m in configuration space by the same amount, which can be described by a 
uniform shift of their indices: M — s- M 4- (rii, 712, 713). 

As far as charge transport is concerned, the important observation is that, since the system remains insulating 
throughout the path, the regions TZm. remain disconnected, so that no charge can flow between them during the 
adiabatic motion. The resulting integrated current flowing through the system during the cycle (which, according to 
Eqs. ^ and ^, measures the change in polarization) can then be inferred from Figs. 2 and 3, which show that the 
regions TZm are labelled by the center of mass of the electrons, and that one can go from one region to the next by 
moving any one electron across the length L of periodic the system. The smallest non-zero change in the average 
polarization along the i-th direction is given by the smallest non-zero shift in the distribution p{Xi) {rii = ±1), and 
is seen to equal \qe\L/V — \qe\/L^, which is the quantum [|l2[|l4| ]. 

It should be noted that strictly speaking the exact quantization of char ge transport for an insulating system 



with periodic boundary conditons only occurs in the thermodynamic limit 25 1. In fact, the quantization was 
established for a finite size only after averaging over all twisted boundary conditions; for a finite size and periodic 
boundary conditions, there are exponentially small corrections ||2^. This is consistent with Kohn's picture that for 
any finite V the regions TZm are not completely disconnected but have an exponentially small overlap, which allows 
for a correspondingly small charge flow between neighboring regions, thus destroying the exact quantization. 

APPENDIX D: THE CASE OF FRACTIONAL FILLING 

The formulas given in the text need to be generalized in order to deal with correlated systems which have in the 
ground state a non-integer number of electrons per primitive cell, and yet are insulating. This can happen in many 
ID models |3^j6^, but it appears that in higher dimensions the Mott transition to an insulating state is usually 
accompanied by a breaking of the symmetry (e.g., a charge-density wave) which restores integer filling [ p3| . Therefore 
we shall restrict the ensuing analysis to ID, as done in Refs. p^ , |36|| . 

Aligia p8[ | has shown that in cases where there is a fractional number of electrons per primitive cell, the limits of 
integration over k in the Berry phase formula for the polarization difference (Eq. [l2| ) given in Ref. need to be 
modified, from which it follows that the "single-point" formulas for the polarization and localization derived in Refs. 
[ p7| , p4t also need to be changed. For the purposes of the present paper it is straightforward to modify the integral over 
k in the definition of InC(a), and from that derive the required modifications to the discrctized numerical formulas. 

Let us consider a ID system with a simulation cell of size L and n/l electrons per unit cell, where n/l is an irreducible 
fraction. Following Ref. ||2^, we modify the cumulant generating function, Eq. ^ as follows: 

InC(a) = 7 / dklnC{k,a). (Dl) 

t .In 



In order to obtain the discretized formulas, we just need to retrace the steps taken in Sections VII A and VII B, with 
the above modification. Discretizing the interval < k < 27r/ /L into a uniform row of J -I- 1 points kj , we find that 
Eq. changes to 
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J-1 



W=^-^Inilnn(<i>^.|<i>fe-,+.). 

7=0 



(D2) 



which agrees with the result quoted in Ref. Similarly, Eq. [79| becomes: 



■In 



7^0 



fc-, I ^fc-y-f 1 ) 



(D3) 



As for the modified "single-point" formulas |^ , it is straightforward to verify, applying the approach of Section VII B 
to the previous two equations, that they are the same as Eqs. |9^ and except for the substitution L L/l. 
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FIGURE 1: Schematic representation of the regions 7?.mi...mjv in the many-electron configuration space where the 
wavefunction of a finite system obeying periodic boundary conditions is locaUzed. The system has hnear dimensions 

a, and the periodic boundary conditions are over a length L >> a, so that the regions 7?,mi...mjv are essentially 
non-overlapping . 



FIGURE 2: Schematic representation of the essentially disconnected regions TZm. in the many-electron configura- 
tion space where the wavefunction of an extended insulator is localized (adapted from Ref. |9|]). For d real-space 
dimensions M — {Mi, ...,Md). Shown is the case of two electrons and d = 1, for which the configuration space is 
{xi,X2)- The system is composed of possibly strongly overlapping units in real space (e.g., a covalent insulator) and 
yet, because it is insulating, in configuration space the wavefunction 5" breaks up into a sum of partial functions ^'mj 
each localized in a region TZm., which have an exponentially small overlap with one another, if the system is large. 



FIGURE 3: Localized distribution p{Xi) along the i-th direction of TV times the electronic center of mass 
(X = X^jLi^j) for a c?-dimensional insulator with N electrons in a periodic volume L"^ (based on Fig. 6 of Ref. 
1^). Although each individual electron coordinate Xj, as well as the electronic charge density (not shown), are spread 
over the whole system, because the system is insulating a localized distribution p{Xi) of width V^iV^i can be uniquely 
defined in terms of the partial wavefunction ^'m in a single disconnected region TZ-m in the dA'^-dimensional configu- 
ration space (Fig. 2). Choosing a different region TZ-^' simply shifts the center of the distribution, if ^ Mi. The 
solid lines correspond to d = 1, for which the peaks coming from different TZ-m do not overlap with one another for 
large L, whereas for d = 3 (dotted lines) they overlap strongly (but the regions TZm are still essentially disconnected, 
for large L.) 



TABLE I. Comparison between the formulas for the Drude weight and for the locahzation length, their relation to the optical 
conductivity, and their asymptotic values in the thermodynamic limit for insulators and conductors at T = 0. 





Drude weight 


Localization length 


Formula in terms of 
twisted boundary conditions 


n _ 1 e'^E(k)\ 


C?(iV) = -^^^|^| 

i 1 CK — 


Relation to conductivity 


Di = — 1 liniij^o ioImaii{uj) 




Asymptotic value 
(A^, ^ oo) 

Insulators 
Non-ideal conductors 
Ideal conductors 


Zero 
Zero 
Finite 


Finite 
Infinite 
Infinite 
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